sigma = 13;
N = 1e5;
y_t = 0:0.1:65;
p_t = y_t/sigma^2;
exp(- y_t.^2 / 2 / sigma^2 );
x_n = rand(1, N);
y_n = sigma*sqrt(-2*log(x_n));
[h, y] = hist(y_n, 30);
inte = sum(h) * (y(2) - y(1));
p = h / inte;
figure(3);
plot(y_t, p_t,'g','LineWidth', 5);
hold on; bar(y, p); 
hold off
xlabel('y'); 
ylabel('p(y)');